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We study the two-dimensional Hubbard model in the weak-coupling regime and compare the 
self-energy obtained from various approximate diagrammatic schemes to the result of diagrammatic 
Monte Carlo simulations, which sum up all weak-coupling diagrams up to a given order. While 
dynamical mean-field theory provides a good approximation for the local part of the self-energy, in¬ 
cluding its frequency dependence, the partial summation of bubble and/or ladder diagrams typically 
yields worse results than second order perturbation theory. Even widely used self-consistent schemes 
such as GW or the fluctuation-exchange approximation (FLEX) are found to be unreliable. Com¬ 
bining the dynamical mean-held self-energy with the nonlocal component of GW in GW+DMFT 
yields improved results for the local self-energy and nonlocal self-energies of the correct order of 
magnitude, but here, too, a more reliable scheme is obtained by restricting the nonlocal contribu¬ 
tion to the second order diagram. FLEX+DMFT is found to give accurate results in the low-density 
regime, but even worse results than FLEX near half-filling. 
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I. INTRODUCTION 

Numerically exact approaches for the solution of corre¬ 
lated lattice models such as the Hubbard model are lim¬ 
ited to one dimension, ’ small lattices,'* weak coupling , 1 
high temperature, or to models with particular sym¬ 
metries and fillings. It is therefore important to de¬ 
velop approximate methods which work in the ther¬ 
modynamic limit, in more than one dimension, and in 
the most interesting range of interactions and densi¬ 
ties. Typically this means interactions comparable to 
the bandwidth and densities close to but not at half 
band filling. One widely used scheme is the dynamical 
mean-field theory (DMFT),' which corresponds to the 
summation of all local self-energy diagrams, via a self- 
consistent impurity construction. This approximation 
becomes exact in the limit of infinite dimensions, ’ as 
well as in the atomic limit and the noninteracting limit. 
It also captures strong-correlation phenomena such as the 
Mott transition. The DMFT approximation however ne¬ 
glects spatial fluctuations and thus cannot be expected to 
capture all the relevant physics in low-dimensional sys¬ 
tems. One possibility is to extend DMFT into a cluster- 
DMFT formalism, which explicitly treats the correla¬ 
tions within some small cluster. Another possibility is to 
implement a diagrammatic expansion around the DMFT 
solution by computing the impurity vertex. “ Both 
approaches are computationally expensive and hence in 
practice are limited to small clusters, or involve the trun¬ 
cation of the diagrammatic expansion to leading order 
terms, or some ladder-type series. Especially in view of 
possible applications to realistic multi-band systems, it 
is thus desirable to devise simpler, computationally less 
demanding schemes. 

One strategy, which has been recently explored in 
simple model contexts, “ is to combine the local 


DMFT self-energy with the nonlocal component of some 
many-body perturbation theory (MBPT) such as second- 
order weak-coupling perturbation theory (lA 2 )) or the 
GW 1 approximation. ’ Alternative schemes, such as 
the combination with the nonlocal self-energy from the 
fluctuation-exchange approximation (FLEX) 2 or the 
T-matrix approximation (TMA ), 2 ' 5 will also be consid¬ 
ered in this work. The advantage of an approach which 
combines DMFT and MBPT at the single-particle (self¬ 
energy) level is that the computational effort is compara¬ 
ble to single-site DMFT and that the extension to multi¬ 
band systems is rather straightforward. The hope is that 
the local self-energy contribution from DMFT captures 
the strong-correlation effects while approximately cor¬ 
rect nonlocal components are introduced by the weak- 
coupling approach. 

In a sufficiently weakly correlated system, the lo¬ 
cal DMFT contribution may not be needed, so that 
self-consistent re-summations of certain classes of weak- 
coupling diagrams, such as bubble and/or ladder dia¬ 
grams, provide an adequate description. While some 
tests of the GW , 2 ’ TMA • or FLEX approaches 
have been published, we still lack a clear picture about 
the importance of the different diagram classes, and the 
beneficial or detrimental effect of self-consistent partial 
re-summations. 

The purpose of this study is to shed some light on these 
issues by benchmarking the approximate self-energies ob¬ 
tained from various MBPT schemes, DMFT and com¬ 
bined MBPT+DMFT approaches against results ob¬ 
tained in diagrammatic Monte Carlo (DiagMC) calcu¬ 
lations, which take into account all diagrams up to a cer¬ 
tain order. More specifically, we focus on the single-band 
Hubbard model on the square lattice 
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with tij = t for i and j nearest-neighbor lattice sites, and 
zero otherwise. The Fourier transform of the hopping 
matrix is hence e*. = — 2t(cos k x + cos k y ). We choose the 
hopping amplitude t = 1 as the unit of energy. Our test 
calculations will be limited to the weak-coupling regime 
U < 4< (half bandwidth), because in this regime con¬ 
verged DiagMC data can be obtained. Such a compari¬ 
son is useful despite this limitation, since a controlled ap¬ 
proximation based on weak-coupling diagrams, or a com¬ 
bination of weak-coupling diagrams and DMFT, should 
behave properly in this limit. 

The paper is organized as follows. In Section II, we 
briefly discuss a number of established approximations 
(DMFT, £( 2 ), GW, TMA, FLEX) and the DiagMC 
method. In Section III, we benchmark the quality of the 
local DMFT self-energy, the local and nonlocal MBPT 
self-energies and various MBPT+DMFT approaches. We 
also study the convergence properties of partial summa¬ 
tions of different classes of weak-coupling diagrams. Sec¬ 
tion IV contains a summary and conclusion. 


II. METHODS 

A. Dynamical mean-field theory 

Dynamical mean-field theory maps a lattice model 
onto a self-consistently defined quantum impurity model 
described by the action 

,1/T 
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where T is the temperature and A(r) is the hybridiza¬ 
tion function. In this approximation the self-energy is 
assumed to be momentum-independent, i.e., E(fe,iw n ) = 
Edmft(*w„). The DMFT self-consistency condition de¬ 
mands that the impurity Green’s function is identical 
to the local lattice Green’s function: f (dfc) G(k,iu n ) = 
Gim P (iu} n ), where /(dfc) denotes a normalized integral 
over the first Brillouin zone. This condition fixes the non¬ 
interacting impurity Green’s function Qq( iui n ), or equiv¬ 
alently the hybridization function A (iw n ) = iui n + /r — 
1/Q 0 (iuj n ), which plays the role of the dynamical mean 
field. In practice, the self-consistent solution is found by 
iterating the following steps (here formulated in terms of 
the “mean field” Go): 

1. Solve impurity model: given G 0 (iuj n ), compute 

G‘ U np (ibj 

n) i 

2. Extract self-energy Edmft/W'AJ = G^iioJn) — 

^imp (® w n)i 

3. DMFT approximation: E (k,iu> n ) = Edmft(*w„), 


4. Compute Gi oc (*w ra ) = J(dk)[iui n + /r — £fc — 
Edmft/A^)] -1 , 

5. DMFT self-consistency: Go 1 fan) = 

Sdmft(*w„) + G loc (iui n ). 

In the present study, the impurity models are solved using 
a strong-coupling continuous-time quantum Monte Carlo 
impurity solver which is numerically exact within sta¬ 
tistical errors. 

The DMFT self-energy corresponds to the sum of 
all one-particle irreducible self-energy diagrams which 
contain only local dressed propagators Gi oc . 59 This 
approximation becomes exact in the limit of infinite 
dimensions. ’ In low-dimensional systems it is a priori 
unclear how important the neglected contributions from 
diagrams with nonlocal propagators are, even for the lo¬ 
cal self-energy. 

B. Weak-coupling approaches 

MBPT encompasses several techniques which, moti¬ 
vated by diagrammatic perturbation expansions, approx¬ 
imate the electron self-energy E at different levels. Meth¬ 
ods like GW or FLEX are frequently considered since 
they can treat spatial fluctuations, are easily imple¬ 
mented and appealing on physical grounds. While the 
truncation of the weak-coupling series for the self-energy 
at the first order yields the Hartree-Fock approximation, 
which for the paramagnetic Hubbard model just amounts 
to a mean-field shift of the chemical potential, the second- 
order approximation, displayed in the top row of Fig. 1, 
includes some nontrivial correlation effects and creates a 
nontrivial frequency and momentum dependence in the 
self-energy. We will denote the second order approxima¬ 
tion with bare propagators by E*' 2 -*. 

Because a systematic computation of the weak- 
coupling expansion to high orders requires rather in¬ 
volved and costly numerical computations (see Sec¬ 
tion IIC), and has only recently become feasible, typical 
approaches to go beyond second order single out specific 
diagram topologies, which may be expected to be domi¬ 
nant in some scenarios, and sum these diagrams analyti¬ 
cally to infinite order by means of Dyson-like equations. 
In addition to the choice of topologies to be included, the 
diagrams can be evaluated with bare propagators Go or 
interacting propagators G. If the self-energy is derived 
from a functional of the self-consistently computed G, 
the approximation can be shown to satisfy certain con¬ 
servation laws. In practice, however, bare expansions 
in terms of Go are often found to be more reliable. ’’ 1 ~ 
(MBPT approaches which involve both bare and inter¬ 
acting propagators have also been proposed, but will 
not be considered here.) 

A further issue is whether the on-site interaction is con¬ 
sidered to act only between different spin species or also 
between identical spins. While the full diagrammatic the¬ 
ory respects the Pauli exclusion principle by construction 
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FIG. 1: (Color online) Illustration of different many-body ap¬ 
proximations to the self-energy. The red wiggly lines represent 
the on-site interaction U. The blue lines with arrows corre¬ 
spond to either bare propagators Go, or (in the case of self- 
consistent perturbation theory) bold propagators G. The first 
(“tadpole”) diagram is the Hartree term. E (2> is the second 
order perturbation theory. Esgw is the spin-dependent GW 
approximation. The spin-independent GW approximation 
in addition contains all the bubble diagrams with an odd num¬ 
ber of interaction lines (not shown). Etma is the T-matrix 
approximation. Finally, Eflex is the fluctuation exchange 
approximations ’ ' 

and is hence oblivious to this choice, diagrammatic ap¬ 
proximations often violate this constraint. For instance, 
the two diagrams in the first row of Fig. 1 constitute all 
terms of the weak-coupling expansion up to third-order 
corrections in U if the interaction only acts between dif¬ 
ferent spins. With a spin-independent interaction, how¬ 
ever, both diagrams would come with a factor of two 
from the spin sum associated with the fermion loops; 
this factor would need to be compensated by a first- and 
a second-order exchange diagram with the same value 
but opposite sign as the shown diagrams. At higher or¬ 
ders many more diagrams need to be included to fully 
compensate terms violating the exclusion principle. We 
therefore adopt the spin-dependent formalism for all the 
approximations shown in Fig. 1. 


In the GW approximation, ’ ’ the self-energy is 
given by the product of the Green’s function G and the 
screened Coulomb interaction W, where only contribu¬ 
tions from the bubble diagrams are considered in the 
calculation of W. In scenarios with long-range Coulomb 
interactions the individual diagrams with bubble inser¬ 
tions are strongly divergent and it is hence essential to 
sum the infinite series into a screened interaction. We 
consider both the self-consistent GW scheme, where all 
propagator lines denote the dressed G, and the “one- 
shot” approach G 0 Wq, where the diagrams are evaluated 
with bare propagators. While most GW calculations as¬ 
sume a spin-independent interaction, this leads to the 
inclusion of W diagrams with an odd number of bare in¬ 
teraction lines, which vanish for the on-site interaction 
of model (1). As this choice (and the neglect of dia¬ 
grams restoring the Pauli principle) effectively removes 
spin-fluctuations, which can be expected to be relevant 
particularly in the vicinity of half filling, we also consider 
the spin-dependent GW approximation (SGW), which 
retains only the even-order diagrams, as illustrated in 
the second row of Fig. 1. The TMA approach, on the 
other hand, sums the series of particle-particle ladder 
diagrams (third row of Fig. 1), which dominate the dia¬ 
grammatic series when the typical inter-particle distance 
is much larger than the range of the interaction. In 
the FLEX approach, finally, bubble, particle-particle 
and particle-hole ladder diagrams are included (bottom 
part of Fig. 1), which means that this approximation 
treats the interaction of electrons via spin, density and 
pairing fluctuations on equal footing. All of these approx¬ 
imations have been widely used to study the properties 
of interacting lattice models or realistic materials in the 
weak-to-intermediate correlation regime. 

The computational steps for the spin-independent GW 
approximation are as follows: 

1. Initialize the self-energy Sci^(fc,iw n ) = 0. 

2. Calculate the Green’s function G(k,ioj n ) = 

1 /K +At- £fc - ?‘Gw(k,iu} n )}. 

3. Calculate the particle-hole polarization function 
Il G w{k,iVm) = ^(T/N k )Y Jq Y, l u lrl G (. c k iu} n)G{q - 
k, iuj n - iv m ). 

4. Calculate the fully screened interaction 
W(k,iv m ) = l/^fc 1 - II Gw{k,iv m )\- For the 
Hubbard model (1), the bare interaction is Vk = U. 

5. Calculate the new self-energy Y>Qw(k,iuj n ) = 
~(T/N k ) Eq G(q,iu n - w m )W{k - q,iv m ). 

6. Go to step 2 until converged results for 
^GwikGwn) and n Gw{k-,iv m ) are obtained. 

Here, uj n denotes a fermionic Matsubara frequency and 
v m a bosonic Matsubara frequency. Nk is the number of 
momentum points in the discretized Brillouin zone. Note 
that in practice we perform the convolutions in the time 
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domain, which allows an efficient treatment of the high- 
frequency components. For the GqWq scheme only one 
pass through steps 1 -5 is performed. 

When the interaction is considered as spin-dependent, 
the equation for W in step 4 should be read as a matrix 
equation in spin-space with a diagonal polarization II and 
an off-diagonal bare interaction v k = Ua x . Its solution 
for the diagonal screened interaction yields 




u 2 u(k, iu m ) 

1 - {UIi{k,iv m )] 2 ' 


( 3 ) 


Additionally, the factor of two in the definition of the 
polarization, coming from the sum over spins, is dropped 
in the spin-dependent case. 

The computational steps for the self-consistent TMA 
calculation are as follows: 

1. Initialize the self-energy StmaC^, = 0. 

2. Calculate the Green’s function G(fc, iuj n ) = 
l/[iu> n + H — e k — S T ma (k, iw„)]. 

3. Calculate the particle-particle polarization function 
n TM A(fe, iv m ) = (T/N k ) J2 q G (<?> iu n )G(k - 
q,w m - iui n ). 


8. Calculate the effective interaction for the 

particle-particle channel V p p {q,iv m ) = 

U/[ 1 + GTIpp(g, iv m )] + U 2 Il pp (q, iv m ). 

9. Calculate the new self-energy Xfle x(k,iu} n ) = 
(T/N k ) E q E Wm [Vph(g, Wm)G(k - q, iu n - iu m ) + 
kp P (*b iVm)G(q k, ii'm i^n)\- 

10. Go to step 2 until Sfle x(k,iui n ) converges. 

Our definitions differ from the original literature - ’ in the 
sign of the particle-hole polarization function II p h, where 
we use the same convention as in the GW scheme, and in 
our inclusion of the Hartree term in the particle-particle 
interaction Vp p , which corresponds to the T-matrix up to 
the correction for the second-order term included in V p h- 
For the non-self-consistent FLEX scheme (FLEXo) only 
one pass through the steps 1-9 is performed. 

Note that in all the above calculations the chemical 
potential fi has to be adjusted self-consistently to en¬ 
sure convergence at the desired density. In our calcula¬ 
tions the ^-summations are discretized in the first Bril- 
louin zone on an equidistant 80 x 80 grid. Furthermore, 
we include the Hartree term in the chemical potential 
rather than the self-energy. In other words we redefine 
the chemical potential and self-energy as 


4. Calculate the T-matrix T{k,iv m ) = —U/[ 1 + 
t/n TM A(fc, Wm)\- 

5. Calculate the new self-energy Etma= 
~{T/N k ) J2 q Yj iVm T{q, Wm)G(q - fc, iv m - iu n ). 

6. Go to step 2 until XTMA(fcAw n ) converges. 

For the non-self-consistent TMA scheme (TMAo) only 
one pass through the steps 1-5 is performed. 

Finally, the procedures for the self-consistent FLEX 
calculation are as follows: 

1. Initialize the self-energy XFLEx(fe, *w n ) = 0. 

2. Calculate the Green’s function G(k,ico n ) = 
l/[iu n + M — £fc — XpLEx(fc, *w„)]. 

3. Calculate the particle-hole polarization function 

n ph (fc, iv m ) = (: T/N k ) Eq G (<r - 

k,iu) n - ivm)- 

4. Calculate the particle-particle polarization function 
n pp (k,iv m ) = {T/N k )J2 q T, ilJJri G(q,iu} n )G(k - 
q,iv m - iu n )• 

5. Calculate the charge susceptibility Xc^AGn) = 

n p h(^ r A i/, m)/[l Gn p i 1 (g, 

6. Calculate the spin susceptibility Xs(?A t 'm) = 

^fph“t“ GHph( q, • 

7. Calculate the effective interaction for 

the particle-hole channel V p \^qG^m) = 

U 2 [§ Xs{q, Wm) + \xMi Wm) + n ph (q, iv m )\. 


=n - Un/2, S' =E - Un/2, (4) 

with n = (rii| + nq.) the number of electrons per site, 
and start all calculations with a “bare” propagator 

G 0 (k , iL 0 n ) =1 /[iu„ + n' - e k \ (5) 

which includes the mean-field effects of the interaction. 
This choice is mostly relevant for one-shot calculations 
and corresponds to the practice in ab initio GW calcu¬ 
lations, which commonly start from a Hartree-Fock or 
density functional solution. ’ 

C. Diagrammatic Monte Carlo 

The DiagMC technique ’ :! 1 evaluates a weak- 
coupling expansion for the self-energy S (k,iui n ) up to 
relatively high orders by means of stochastic sampling. 
In contrast to the approximate schemes discussed above, 
all diagram topologies are included. A few examples of 
diagrams neglected in the previous schemes are shown 
in Fig. 2. While at least FLEX includes all topologies 
occurring up to third order, the majority of fourth-order 
diagrams is already neglected. For higher orders, only an 
exponentially small fraction of the diagrams at a given 
order is included in approximate methods such as GW, 
TMA or FLEX. 

Both the sums over diagram orders and topologies, and 
the integrals over internal variables are sampled using a 
Monte Carlo procedure. By restricting the sampling pro¬ 
cess to one-particle irreducible diagrams the self-energy is 
computed directly and can then be inserted into Dyson’s 
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FIG. 2: (Color online) Examples of diagram topologies not 
contained in any of the expansions shown in Fig. 1. (a) Self¬ 
energy insertions on internal propagator lines. [These are ac¬ 
counted for in self-consistent schemes which use an expansion 
in terms of the interacting propagator.] (b) Ladders with 
crossed rungs, (c) Topologies with more complex vertex cor¬ 
rections. 


FIG. 3: (Color online) Convergence of the weak-coupling 
series for the local self-energy with diagram order TV* and 
comparison to the DMFT self-energy (solid horizontal lines). 
Shown are the real (black) and imaginary (blue) parts at the 
lowest Matsubara frequency loo = 17tT for two different fill¬ 
ings h = 0.4, 0.8 and two values of the interaction strength 
U = 2, 4. The temperature is T = 0.1 in both cases. 


equation to obtain a single-particle propagator G{k, ico n ) 
corresponding to an infinite number of diagrams, com¬ 
posed of arbitrary combinations of the explicitly sam¬ 
pled self-energy diagrams. The only systematic error 
consists in a cutoff of the diagrammatic series at order 
A r *, i.e., the weak-coupling diagrams are generated for or¬ 
ders N < IV*. Such a cutoff must be introduced because 
the average sign in the Monte Carlo sampling vanishes 
exponentially with growing diagram order. By varying 
IV* and monitoring the convergence of the self-energy, 
the accessible parameter regime can be determined and 
the errors can be controlled. We use an expansion in 
terms of bare propagators which is typically found to con¬ 
verge towards the correct solution in the weak-coupling 
regime U < 4f, wherever numerically exact benchmarks 
are available. ’ 11 


III. RESULTS 

In the following we compare the self-energies ob¬ 
tained from DMFT, several weak-coupling approxima¬ 
tions and MBPT+DMFT schemes to the accurate and 
well-controlled DiagMC self-energy. We concentrate on 
the nontrivial part of the self-energy, E', obtained after 
subtraction of the Hartree contribution [see Eq. (4)]. 


A. Local self-energy 

We start by benchmarking the local self-energy ob¬ 
tained within DMFT. Figure 3 plots the lowest Mat¬ 
subara frequency component of the local self-energy 
EJoc(iwo) = /(d/c) E'(fc, iujo) calculated by the Di¬ 
agMC method up to order IV* = 7 and compares it 
to Edmft(*^o) for two different site occupancies n = 
0.4 and 0.8 and interaction strengths U = 2 and 4. 
We find that both the real and imaginary parts are 
quite accurately reproduced by DMFT: ImEDMFT(*wo) 
agrees with the DiagMC result within error bars, while 
ReEDMFT(*wo) deviates by less than 10%. 

While the momentum-dependence is neglected, DMFT 
can capture a nontrivial frequency dependence of the self¬ 
energy. Figure 4 shows the comparison of this frequency 
dependence to the DiagMC results for the same param¬ 
eter sets. We see that DMFT also predicts the correct 
frequency dependence of the local self-energy, with max¬ 
imum relative deviations of less than 10%. 

We next consider the local component of the self¬ 
energy obtained from different weak-coupling approxima¬ 
tions. Figure 5 shows the comparison of the E^ 2 ), GW, 
S GW, TMA and FLEX results to DiagMC, for the same 
parameters U = 2, 4 and n = 0.4, 0.8. While none of the 
weak-coupling approximations are as accurate as DMFT, 
the E^ 2 ) (and to a lesser extent the SGW) approximation 
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n =0.8 




FIG. 4: (Color online) Frequency dependence of the local 
self-energy from DiagMC (black and blue dots for the real 
and imaginary parts, respectively) and DMFT (solid lines) 
for the same systems as in Fig. 3. DiagMC error bars cover 
the results with the three largest cut-off orders N, = 5, 6, 7. 


reproduces the exact results rather well. FLEX gives rea¬ 
sonable estimates for the real part but can significantly 
overestimate the imaginary part, especially near half fill¬ 
ing. GW and the TMA yield poor estimates of either the 
real or imaginary part. Not surprisingly, the quality of 
the TMA decreases with increasing interaction strength 
and away from the dilute limit. The GW approxima¬ 
tion, on the other hand, tends to strongly overestimate 
the self-energy for weak interactions. Based on these re¬ 
sults, we must conclude that schemes involving partial 
summations of diagrams are less reliable than the sim¬ 
ple approximation. Additionally, X^ 2 ) is the only 
weak-coupling approximation that correctly captures the 
exact asymptotic behavior of the self-energy at high fre¬ 
quencies, X( oc (»« n ) = U 2 n a (l-fi a )/(iuj n )+0(l/(iu} n ) 2 ), 
whereas all the other schemes significantly over- or un¬ 
derestimate the coefficient of this tail (not shown). 

The poor performance of the spin-independent GW 
approximation may at first seem surprising given its 
widespread and successful use in electronic structure cal¬ 
culations. However, as discussed above, all the other 
schemes assume a spin-dependent interaction and con¬ 
tain the correct first- and second-order terms of the weak- 
coupling expansion for a local interaction, whereas GW 
overestimates the latter term by a factor of two. In ab 
initio calculations of weakly correlated materials, where 
GW is primarily used, local interactions, and hence vi¬ 
olations of the exclusion principle, may be expected to 
be less relevant. With a non-local interaction, also the 
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FIG. 5: (Color online) Frequency dependence of the local self¬ 
energy in the X (2 1, GW, SGW, TMA and FLEX approxima¬ 
tion compared to the same DiagMC results as shown in Fig. 4. 
The top four panels plot the real parts, and the bottom four 
panels the imaginary parts of the self-energy. 


other schemes would need to be supplemented by addi¬ 
tional exchange diagrams in order to correctly capture 
all weak-coupling contributions. 

Here and in the following we concentrate on the self- 
consistent versions of MBPT - except for the simple X*- 2 -* 
approximation for which we show the one-shot result. 
In the parameter regime considered here, the difference 
between one-shot and self-consistent calculations is small 
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for n = 0.4, while there can be significant differences for 
n = 0.8. An explicit comparison between one-shot and 
self-consistent results for the data of Fig. 5 is shown in 
Appendix A. 

In view of these results, the idea of replacing the local 
component of the MBPT self-energy by the more reli¬ 
able DMFT self-energy appears to be reasonable. But 
before we investigate how this replacement affects differ¬ 
ent self-consistent schemes, we take a look at the nonlocal 
components of the self-energy. 


B. Nonlocal self-energy 

Since in the weak-coupling regime considered here the 
nonlocal self-energy £^ decays rapidly with the distance 
* — j\, and it is computationally expensive to obtain Di¬ 
agMC data with small error bars, we restrict the tests of 
the nonlocal components to the nearest-neighbor contri¬ 
bution £„„ and the next-nearest-neighbor contribution 
£„„„. Figure 6 shows the frequency dependence of these 
components, again for U = 2, 4 and n = 0.4, 0.8. Both 
the real (black) and imaginary (blue) parts are plotted 
in the same panel, and the error estimates of the Di- 
agMC results are indicated by grey and blue shading. 
We have estimated the systematic uncertainty on the Di- 
agMC data by considering the results for the four largest 
cutoffs, while the stochastic uncertainty is estimated from 
64 independent runs. 

By comparing the y -axis scales in Fig. 6 to the corre¬ 
sponding plots for the local component of the self-energy 
(Fig. 5) we see that the £„ n and £ raTm are at least a fac¬ 
tor of ten smaller. While the weak-coupling approxima¬ 
tions produce nonlocal components of the correct order 
of magnitude, the relative errors are large. None of the 
weak-coupling approximations gives reliable results for 
both the real and imaginary parts. While FLEX seems 
to work well for the imaginary part of £ ran , it gives poor 
results for the real part and for T, nnn . GW and the TMA 
do not produce very inaccurate results but they are not 
systematically better than £( 2 ^. To avoid overcrowding 
the figure, we have not plotted the SGW results, which 
are typically between those of GW and FLEX. As for 
the local self-energy, we conclude that there seem to be 
no obvious benefits from partially summing diagrams be¬ 
yond the second order. 


C. Combinations of DMFT with weak-coupling 
approximations 

Since the DMFT approximation provides a good de¬ 
scription for the dominant local part of the self-energy, 
and weak-coupling perturbation theories produce at least 
a reasonable estimate of the nonlocal components, it 
is tempting to combine the two approaches. Indeed, 
such methods have been proposed many years ago, in 
particular the combination of and DMFT and 
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FIG. 6: (Color online) Frequency dependence of the nonlo¬ 
cal self-energy for nearest- and next-nearest-neighbor sites in 
the £ (2 \ GW, TMA and FLEX approximation compared to 
DiagMC results. The shaded bands cover stochastic lcr er¬ 
rors around the DiagMC results with the four largest cutoff 
orders A* = 4,..., 7. The top four panels plot the self-energy 
for nearest neighbors A r = (1,0) and the bottom four panels 
for next-nearest neighbors along the diagonal of the square 
lattice A r = (1,1). 


the combination of GW and DMFT. These methods 
have been designed in particular to treat models with 
long-ranged Coulomb interactions, based on an extended 
DMFT (EDMFT) formalism, ’ ’ ’ and because of re- 


























cent methodological advances related to impurity prob¬ 
lems with dynamically screened interactions, 45,11 there 
has been a revival of interest in these approaches. ' 6> 
The same techniques can also be applied to model (1) 
with only an on-site Hubbard repulsion. We will consider 
here the S( 2 )+DMFT, GW+DMFT and FLEX+DMFT 
schemes, in which the lattice self-energy is approximated 
as 


^MBPT+DMFT 

-‘jk 


(*w„) (io, n )J jfc 


, vMbpt 
+ ^jk 


(w„)( l-d jk ). (6) 


We have also implemented TMA+DMFT, but will not 
show these results, because they do not change the main 
conclusions. Note that there are various ways of pre¬ 
venting the double-counting of diagrams. Equation (6) 
corresponds to the simplest approach, i.e., the removal of 
all the local MBPT self-energy diagrams. This double¬ 
counting scheme also removes diagrams with nonlocal 
propagators, which are not included in the DMFT self¬ 
energy. An alternative way of combining the DMFT and 
MBPT diagrams is 


yiMBPT+DMFT / • x _ v 
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where E^ 4BPT [G J y](jw r i) denotes the subset of E^ 4BPT di¬ 
agrams which contains only local propagators G l: j . We 
have tested both double counting schemes, but for the pa¬ 
rameter sets considered, the differences are rather small. 
We will show the results for the self-energy (6) , and com¬ 
ment in the text on the effect of the alternative scheme 
(7), where appropriate. 

Because the MBPT+DMFT calculations are done self- 
consistently, it is not easy to identify the subsets of dia¬ 
grams summed up by these schemes. However, as can be 
seen in Fig. 7, the local X in the GW + DMFT approx¬ 
imation reproduces the DiagMC result very well. The 
imaginary part agrees with DiagMC within error bars, 
and is thus even more accurate than the DMFT result 
(Fig. 4), while the accuracy of the real part is compara¬ 
ble to DMFT. Since the real part is very sensitive to the 
value of the chemical potential, some of these differences 
may be explained by the uncertainty in the self-consistent 
calculation of p,. 

In Refs. 16 and 17 it was found that the combina¬ 
tion of GW and EDMFT makes the system more corre¬ 
lated, compared to EDMFT. This conclusion was based 
on an extended Hubbard model calculation at half filling, 
with U = 8t and nearest neighbor Coulomb repulsion 
V > 0.8 1. Comparing Figs. 7 and 4 we find the oppo¬ 
site effect in the simple Hubbard model away from half 
filling: the imaginary part of the self-energy is slightly 
reduced by adding the nonlocal GW self-energy, which 
means that the system becomes less correlated. This dif¬ 
ference may be due to the fact that we consider here a 
less correlated system, a system away from half filling, 
or it may indicate that the enhanced correlations in the 



FIG. 7: (Color online) Frequency dependence of the local self¬ 
energy obtained by £®+DMFT (solid lines), G1F+DMFT 
(dashed lines), and FLEX+DMFT (dotted lines) compared 
to the same DiagMC data as in Fig. 4. 


previous GLF+EDMFT studies result from a nontrivial 
interplay between the nonlocal self-energy and nonlocal 
screening. In any event, it seems that the addition of 
the nonlocal GW self-energy can both increase or de¬ 
crease the local correlations, depending on the parameter 
regime. 

FLEX+DMFT gives improved local self-energies com¬ 
pared to DMFT for U = 2, and for [7 = 4, n = 0.4, 
but the result for U = 4, n = 0.8 is significantly less 
accurate than the DMFT prediction. (With the alter¬ 
native double counting scheme (7), the real part of the 
self-energy is improved at low Matsubara frequencies, but 
the imaginary part is overestimated.) Apparently, close 
to half filling, the feedback from the inaccurate nonlocal 
FLEX self-energy has a detrimental effect on the local 
self-energy. 

While the differences to GW +DMFT are not very sig¬ 
nificant, the simple E^+DMFT scheme yields the most 
accurate estimates of the local self-energy, for both inter¬ 
actions and fillings. 

Figure 8 compares the nonlocal self-energy compo¬ 
nent E nn obtained from the E^+DMFT, GW+DMFT 
and FLEX+DMFT calculations to the DiagMC re¬ 
sults. The comparison between the MBPT results and 
MBPT+DMFT are shown in Appendix B. These re¬ 
sults illustrate how the self-consistent feedback of the 
DMFT self-energy into the MBPT scheme affects the 
nonlocal self-energy. In the case of E^+DMFT and 
GLF+DMFT, the change with respect to the nonlocal 
El 2 ) and GW self-energy is small and there is no sys- 
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FIG. 8: (Color online) Frequency dependence of the nearest- 
neighbor self-energy obtained by the £ <2 '+DMFT (solid 
lines), GW+DMFT (dashed lines) and FLEX+DMFT (dot¬ 
ted lines) schemes compared to the same DiagMC data as in 
the upper half of Fig. 6. 


tematic improvement of the nonlocal components. For 
FLEX+DMFT, the conclusion is similar in the case of 
U = 2 and U = 4, n = 0.4, while for U = 4, n = 0.8 
FLEX+DMFT is significantly less accurate than FLEX. 
(While the double counting scheme (7) improves the re¬ 
sults somewhat, both the real and imaginary parts of 
X nn are still significantly overestimated.) Hence, in 
the parameter regime where MBPT is not too inac¬ 
curate, the local self-energy is apparently improved in 
the MBPT+DMFT approach, while the nonlocal com¬ 
ponents of X are almost unchanged, and do not system¬ 
atically benefit from the additional local self-energy di¬ 
agrams in the nonlocal propagators. If the MBPT re¬ 
sult deviates strongly from the correct solution, as is the 
case with FLEX in the intermediate coupling regime close 
to half filling, then the self-consistent combination with 
DMFT has detrimental effects on both the local and non¬ 
local components of the self-energy. 


D. Relevant diagrams 

As discussed in Section IIB, a basic assumption un¬ 
derlying approximate schemes such as GW and FLEX 
is that specific diagram topologies with a rather simple 
structure contain the relevant physics, at least in cer¬ 
tain parameter regimes, such that the summation can be 
restricted to a tractable subset. In order to test this as¬ 
sumption and possibly identify the relevant subsets, we 


have implemented a classification scheme for the sam¬ 
pled diagrams in our DiagMC code. This allows us to 
check, order by order, the respective contributions from 
GW-type bubble diagrams or the particle-particle (pp) 
and particle-hole (ph) ladders included in the TMA and 
FLEX approximations. In addition, we consider the class 
of generalized ladder diagrams (“X-ladders” for brevity) 
that includes not only the pp- and ph-ladders but also 
those with crossed rungs, some examples of which are 
displayed in Fig. 2 (b). 

Here, we concentrate on the case U = 4 and study 
the evolution of Xi oc (muo) with increasing diagram or¬ 
der. We first focus on the bare expansion in terms of 
the noninteracting propagator Gq. The left panels of 
Fig. 9 show data for n = 0.4, with the solid black curve 
corresponding to the DiagMC result which sums up all 
diagram topologies. The other curves correspond to the 
above-mentioned families of diagrams and their combina¬ 
tions. We note that the “bubble” diagrams correspond to 
those included in a spin-dependent Go Wo calculation and 
the “pp-ladder” to a one-shot TMAo scheme, while the 
“bubble+ladders” curves contain exactly the topologies 
included in a one-shot FLEX calculation. We indicate 
the results of these one-shot calculations (with the same 
chemical potential as used in the corresponding DiagMC 
simulation) with colored arrows on the y-axis. 

We see that both the bare particle-particle and 
particle-hole ladders start to deviate significantly from 
the exact result for orders > 3, albeit in opposite ways. 
The bare bubble series seems to be slightly better be¬ 
haved although it tends to worsen rather than improve 
the second-order result, in agreement with the findings 
for the SGoWo approximation. While the combination 
of the particle-particle and particle-hole ladders does not 
help much, the inclusion of diagrams with crossed rungs 
in X-ladders does improve the result. This finding is 
consistent with the intuition of Bickers and White, 1 who 
suggested that ladders with crossed rungs should strongly 
renormalize the particle-particle and particle-hole ladder 
contributions, and argued that one should therefore work 
with a renormalized U. (It should be kept in mind that 
the X-ladder class of diagrams cannot be summed ana¬ 
lytically via a Dyson equation.) At least for n = 0.4, the 
sum of bubbles and X-ladders yields a self-energy which 
is relatively close to the exact result for the diagram or¬ 
ders considered here. 

The situation gets worse closer to half filling (n = 0.8, 
see middle panels in Fig. 9). Here, the “bubble+X- 
ladders” result deviates strongly from the full series, at 
least for the real part of the self-energy. Also the other 
diagram families either converge to wrong values or show 
no sign of convergence up to the seventh order. This in¬ 
stability is also evident in the FLEX calculations, which 
need to be initialized with a chemical potential corre¬ 
sponding to a lower filling in order to avoid diverging sus¬ 
ceptibilities in the first iteration. Consequently, there are 
no FLEX 0 results indicated in the central panels. Over¬ 
all, it is clear that none of these families of diagrams 
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[7=4, n =0.8, bare diagrams (7=4, n=0.8, skeleton diagrams 



FIG. 9: (Color online) Convergence of the local self-energy with diagram order for the full series, sampled in DiagMC, and 
various subclasses of diagrams. See main text for an explanation of the different diagram classes. The upper (lower) row shows 
the real (imaginary) part at the lowest Matsubara frequency iui o = inT. The left and center columns correspond to the bare 
series at two different densities, while the right column shows the skeleton series for the same parameters as the central panels. 
The black arrows in the right panel show the converged DiagMC results, as estimated from the bare series. 


yields a systematically better approximation of the lo¬ 
cal self-energy than the second-order contribution. 
Apparently, the cancellation effects among higher order 
contributions are so subtle that essentially all diagram 
topologies must be considered, and the restriction to a 
subset of ladder or bubble type diagrams cannot be jus¬ 
tified. This is further corroborated by the observation 
that all the subclasses converge, if at all, far less regularly 
at large orders than the sum of all topologies. Even the 
X-ladders class, which grows exponentially with diagram 
order, exhibits seemingly erratic kinks beyond the fifth 
order, which are apparently canceled by other diagrams, 
since they are not visible in the sum of all diagrams. 

One may wonder whether the situation can be im¬ 
proved by considering only two-particle irreducible skele¬ 
ton diagrams and replacing the bare propagators Go by 
self-consistently computed interacting Green’s functions 
G. In order to check this hypothesis we conducted a 
DiagMC sampling of skeleton diagrams where the prop¬ 
agators are dressed with the self-energy obtained from a 
previous sampling of the bare series up to sixth order. 
While such self-consistent calculations sum up more dia¬ 


grams, the right panels of Fig. 9 show that the boldified 
diagrammatic series converges more slowly than the bare 
series. This observation is in accord with the recent re¬ 
sults of Ref. 33, where the skeleton series was found to 
converge very slowly at intermediate interaction U ~ At. 
At larger interaction U At. the skeleton series is even 
reported to converge to an unphysical solution, whereas 
the bare expansion shows no such pathological behavior. 
For the shown parameters the bold X-ladders result is 
close to DiagMC, but this good agreement appears to be 
accidental since the corresponding curves at other fre¬ 
quencies significantly deviate from each other, with the 
X-ladders seemingly converging to incorrect values (not 
shown). 


IV. CONCLUSIONS 

We have performed a systematic study of the accuracy 
of various approximate diagrammatic schemes for the so¬ 
lution of the 2D Hubbard model. By comparing the self¬ 
energies obtained from widely used MBPT approaches 
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and DMFT to the well-controlled DiagMC result we were 
able to assess the quality of the approximations in the 
weak-coupling regime. We have also measured order by 
order the contribution of different diagram classes in or¬ 
der to track their convergence properties. The main con¬ 
clusion is that none of the conventional schemes like GW, 
TMA or FLEX, which sum up bubble and/or ladder dia¬ 
grams, provides a systematic improvement over the sim¬ 
ple approximation, and in fact often yield consid¬ 
erably less accurate results. The systematic bias and/or 
the erratic convergence properties of these schemes with 
diagram order indicate that the corresponding small sub¬ 
classes of diagrams do not capture the dominant contri¬ 
butions to the self-energy, and that the corrections from 
the neglected diagrams are significant. Even by consider¬ 
ing additional diagram topologies such as X-ladders, we 
were not able to identify a ‘relevant subset’ of diagrams. 
It thus appears that in general, the partial summation 
of ladder or bubble type diagrams is not a valid approx¬ 
imation, because essentially all diagram topologies are 
relevant. At least in the weak-coupling regime, stopping 
at the second order (X^ 2 )) is more reliable than perform¬ 
ing uncontrolled summations. Although we cannot ac¬ 
cess the intermediate and strong-coupling regime with 
DiagMC, it seems unlikely that a weak-coupling based 
MBPT approach which is found to be unreliable in the 
weak-coupling regime can be trusted in the more strongly 
correlated regime. While MBPT methods have been em¬ 
ployed by many groups to study transition metal and 
actinide compounds, ’ “ our findings put a question 
mark behind the use of GW or FLEX (both the one-shot 
and self-consistent variants) in studies of lattice models 
or materials with substantial local correlations. 

For the local part of the self-energy, the DMFT ap¬ 
proximation, which is nonperturbative and sums all dia¬ 
grams made from local propagators, provides a good ap¬ 
proximation. This class of diagrams can however not be 
summed by a simple Dyson-type equation, but requires 
a self-consistent impurity model calculation. At least 
in the weak-coupling regime, where the nonlocal compo¬ 
nents of the self-energy are small, and as we have shown 
are reasonably described by many-body perturbation ap¬ 
proaches such as X^ 2 ) or GW , it makes sense to combine 
the two approaches by adding the nonlocal component of, 
e.g., the GW self-energy to the local DMFT self-energy. 
We have tested several MBPT+DMFT schemes and 
found that for X( 2 )+DMFT and GW+DMFT the feed¬ 
back from the nonlocal component in the self-consistency 
loop improves in particular the local self-energy, which 
becomes very accurate. The nonlocal components are not 
systematically improved compared to the pure MBPT re¬ 
sult, but of comparable accuracy. In FLEX+DMFT, the 
inaccuracy of the FLEX contribution near half filling can 
lead to self-energies which are significantly less accurate 
than the DMFT prediction. 

While GfF+DMFT has been found to underestimate 
the /c-dependence of the self-energy in the intermedi¬ 
ate coupling regime, 1 '’’ 1 compared to cluster DMFT 


calculations, ;,J ’ this result is not really surprising. The 
GW method has been primarily designed to capture 
the effect of screening from long-ranged Coulomb inter¬ 
actions and to avoid divergences in the diagrammatic 
expansion in terms of the bare Coulomb interaction. 
This is very important for the proper description of 
materials, : ’ but does not play a role in the Hub¬ 
bard model with purely on-site interactions considered 
in this study. The main target for GW+DMFT and re¬ 
lated approaches is thus the realistic simulation of (three- 
dimensional) compounds, where the fc-dependence can be 
expected to be small, while the effect of dynamical screen¬ 
ing may be significant. In this case, the DMFT loop 
should also involve a self-consistent calculation of the 
screened interaction ’ and the impurity model must 
be extended to one with retarded density-density inter¬ 
actions, or even retarded spin-flip terms. 

In this work, we have focused on methods which 
combine DMFT and MBPT at the single-particle (self¬ 
energy) level. These should be distinguished from more 
sophisticated, but also computationally much more ex¬ 
pensive methods which attempt such a combination at 
the two-particle (vertex) level. Since the Coulomb in¬ 
teraction is a two-particle operator, the latter approach 
may be expected to yield better results and access to 
the particularly interesting intermediate coupling regime. 
Systematic tests of methods such as the dynamical ver¬ 
tex approximation or dual fermion based schemes against 
DiagMC results could provide valuable insights into the 
virtues and limitations of these vertex based approaches. 
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Appendix A: Self-consistent vs. one-shot calculations 

In the parameter regime considered, the differences be¬ 
tween one-shot and self-consistent weak-coupling approx¬ 
imations are typically rather small compared to the dif¬ 
ferences between different diagrammatic approximations. 
Figure 10 shows the results for the local self-energy com¬ 
puted using the X/ 2 ), TMA, GW, SGW and FLEX ap¬ 
proximations for U = 2, 4 and h = 0.4, 0.8. For reference, 
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FIG. 10: (Color online) Comparison of self-consistent (broken lines) and one-shot (solid lines) results for the real and imaginary 
parts of the local self-energy from the various weak-coupling approximations. 


the respective DiagMC result is indicated by a grey band 
in each panel. For h = 0.4, the self-consistent resumma¬ 
tion changes the self-energy only slightly. For n = 0.8 
and weak interaction (U = 2) the differences between 
the one-shot and self-consistent schemes are smaller than 
those between the various weak-coupling approximations. 
In the vicinity of half filling and for stronger interaction 
(U = 4, n = 0.8 in Fig. 10), the difference becomes sig¬ 
nificant. Here in particular the spin-dependent GW and 
FLEX approximations are close to a pole in the expres¬ 
sion for the effective interaction and therefore very sen¬ 
sitive to changes in the polarization. Dressing the prop¬ 
agator reduces the polarization’s magnitude and hence 
moves the expression away from the pole, reducing the 
resulting self-energy. 

Appendix B: Effect of DMFT corrections on the 
nonlocal self-energy 

Figure 11 compares the nonlocal self-energy obtained 
from GW, and FLEX to the corresponding results 


produced by the combinations of these MBPT methods 
with DMFT. We see that the inclusion of additional local 
diagrams in the MBPT propagators has only a moderate 
effect - again with the exception of FLEX in the vicin¬ 
ity of half filling and for stronger interactions - and does 
not systematically improve the result. While the double 
counting scheme (7) improves the nonlocal self-energy for 
FLEX+DMFT, the results for both the real and imag¬ 
inary parts are still significantly larger than for FLEX, 
and hence less accurate. 
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n =0.4 
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n =0.8 





FIG. 11: (Color online) Comparison of the self-energy for nearest neighbor sites from E* 2 -* (top left), GW (top right) and 
FLEX (bottom) with the results from the corresponding MBPT+DMFT approaches. In the FLEX case, we have multiplied 
the results for the lower filling n = 0.4 by ten, to make them visible on the same scale as those for n = 0.8. 
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